## Step 4: Infer rho/eta consistent with the distribution of growth rates -----

# Read in CE rates from NP2003 (RW, MW), Groom et al. 2007 (SS), 
# and Freeman et al 2015 (AADL)
ce_rates = array(NA, dim=c(500, 6, 5), 
                 dimnames=list(1:500, c('BR','LVL','RW','MR','SS','AADL'),
                               c('1.5%','2%','2.5%','3%','5%')))

load(data.dir%&%'/BR_term_structures.RData')
ce.BR = ce.BR[-1,-c(1,6)] # drop year 1 which is the near-term rate by construction
sum(ce.BR)
dimnames(ce.BR)
# ce.BR = ce.BR[1:tt,,]
dim(ce.BR)

ce_rates[1:nrow(ce.BR),'BR',] = ce.BR

ce_rates[,'RW','1.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate1pt5pct_RW.csv')[,1]
ce_rates[,'RW','2%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pct_RW.csv')[,1]
ce_rates[,'RW','2.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pt5pct_RW.csv')[,1]
ce_rates[,'RW','3%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate3pct_RW.csv')[,1]
ce_rates[,'RW','5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate5pct_RW.csv')[,1]

ce_rates[,'MR','1.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate1pt5pct_MR.csv')[,1]
ce_rates[,'MR','2%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pct_MR.csv')[,1]
ce_rates[,'MR','2.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pt5pct_MR.csv')[,1]
ce_rates[,'MR','3%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate3pct_MR.csv')[,1]
ce_rates[,'MR','5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate5pct_MR.csv')[,1]

ce_rates[,'LVL','1.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate1pt5pct_LVL.csv')[,1]
ce_rates[,'LVL','2%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pct_LVL.csv')[,1]
ce_rates[,'LVL','2.5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate2pt5pct_LVL.csv')[,1]
ce_rates[,'LVL','3%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate3pct_LVL.csv')[,1]
ce_rates[,'LVL','5%'] = read.csv(data.dir%&%'/Term structure data/NP2003/saveCErate5pct_LVL.csv')[,1]

class(ce_rates[,'SS',])

# Note: groom et al. and freeman et al. rates are forward rates.
# Convert to average rates.
avg.to.fwd = function(r) {
  # Takes an average rate and derives the forward rates.
  # Formula is fwd[t] = avg[t]*t - avg[t-1]*(t-1), with avg[0] = 0
  # avg[1] = fwd[1]
  # avg[2] = (fwd[1]+fwd[2])/2
  # ...
  # avg[t-1] = cumsum_1^t-1(fwd[i]) / (t-1)
  # avg[t] =   cumsum_1^t  (fwd[i]) / t
  # So:
  # fwd[t] = avg[t]*t - avg[t-1]*(t-1)
  # where we define avg[0] = 0, so fwd[1] = avg[1]*1
  tt = 0:length(r)
  r0 = 0 # avg[1] = 0
  return(diff(c(r0, r)*tt))
} 

fwd.to.avg = function(r) {
  result = cumsum(r)/(1:length(r))
  return(result)
} 

# Check that f(finverse(r))==r
fwd.to.avg(avg.to.fwd(ce_rates[,'RW','3%'])) / ce_rates[,'RW','3%']

library(openxlsx)
groom.et.al = read.xlsx(data.dir%&%'/Term structure data/Discount rates from Groom et al.xlsx', sheet=2)[,-1]
names(groom.et.al) #= round(as.numeric(names(groom.et.al)), 3)
groom.et.al = groom.et.al[,-which(names(groom.et.al)=="0.04")] # not the 4% column (col#4)
head(groom.et.al)
groom.et.al = apply(groom.et.al, 2, fwd.to.avg)

ce_rates[1:nrow(groom.et.al),'SS',] = as.matrix(groom.et.al) 

freeman.et.al = read.xlsx(data.dir%&%'/Term structure data/Discount rates from Groom et al.xlsx', sheet=1)[,-1]
names(freeman.et.al)
freeman.et.al = freeman.et.al[,-which(names(freeman.et.al)=="AADL(4%)")] # not the 4% column (col#4)
head(freeman.et.al)
freeman.et.al = apply(freeman.et.al, 2, fwd.to.avg)

ce_rates[1:nrow(freeman.et.al),'AADL',] = as.matrix(freeman.et.al) 

tt = nrow(global_growth_percap)
tt # 280 observations in MWS. 
# Keep only the first 280 of the ce rates
ce_rates = ce_rates[1:tt,,]
dimnames(ce_rates)
